////////////////////////////////////////////////////////////////////////
//irmm_Raster.h
//This is the description of classes of image library including
//======================================================================
//		Gaussian filtering
//		segmentation of color image using mean shift algorithm
//=======================================================================
//	Author: Zhuowen Tu
////////////////////////////////////////////////////////////////////////
#ifndef _IRMM_RASTERCOLOR_H
#define _IRMM_RASTERCOLOR_H

#include "mcRaster.h"
#include "Vector.h"

// List of color components; the possible cases are:
#define RedColor		(1<<0)
#define GreenColor		(1<<1)
#define BlueColor		(1<<2)
#define Lightness		(1<<3)
#define Ustar   		(1<<4)
#define Vstar   		(1<<5)
#define FirstChannel    (1<<6)
#define SecondChannel   (1<<7)
#define ThirdChannel    (1<<8)

#define LUV2RGB         10
#define YCrCb2RGB       11


template<class Tp>
class ColorRGB
{
public:
	Tp r;
    Tp g;
    Tp b;

	ColorRGB(){r=0;g=0;b=0;};
	ColorRGB(Tp ir,Tp ig,Tp ib){r=ir;g=ig;b=ib;};

	inline void Set(Tp rc,Tp gc,Tp bc){r=rc;g=gc;b=bc;};
};

class ColorLUV
{
public:
	int	L;
    int	U;
    int	V;

	ColorLUV(){L=0;U=0;V=0;};
	ColorLUV(const int iL, const int iU, const int iV){L=iL; U=iU; V=iV;};

	inline void Set(const int iL, const int iU, const int iV){L=iL; U=iU; V=iV;};
};

#define MAXV 256

////////////////////////////////////////////////////////////
//this is the class of color image with three channels in r, g, b or
//L, U, V spaces
template<class Tp>
class Raster_Color
{
public:
	// construction
	Raster_Color();
	~Raster_Color();

	// operators
	Raster_Color& operator=(Raster_Color<uchar> &I);
	Raster_Color& operator=(Raster_Color<double> &I);
	void Copy(Raster_Color<uchar> &I, irmm_Rect& rect);

	void	Free_histogram_space(void);

	void SetDimension(const int nrows, const int ncols);

	virtual void	LoadPpmPGm(const char *pstr_name);
	virtual bool	SavePpmPGm(const char *pstr_name, int iff);

	// access
	inline Tp&		R(const int Row, const int Col);
	inline Tp&		G(const int Row, const int Col);
	inline Tp&		B(const int Row, const int Col);
	inline Tp&		R(const long index);
	inline Tp&		G(const long index);
	inline Tp&		B(const long index);

	inline int&		L(const int Row, const int Col){return LUV[0](Row, Col);};
	inline int&		U(const int Row, const int Col){return LUV[1](Row, Col);};
	inline int&		V(const int Row, const int Col){return LUV[2](Row, Col);};
	inline int&		L(const long index){return LUV[0](index);};
	inline int&		U(const long index){return LUV[1](index);};
	inline int&		V(const long index){return LUV[2](index);};

	//add by wtf
	inline int&		Y(const int Row, const int Col){return YCrCb[0](Row, Col);};
	inline int&		Cr(const int Row, const int Col){return YCrCb[1](Row, Col);};
	inline int&		Cb(const int Row, const int Col){return YCrCb[2](Row, Col);};
	inline int&		Y(const long index){return YCrCb[0](index);};
	inline int&		Cr(const long index){return YCrCb[1](index);};
	inline int&		Cb(const long index){return YCrCb[2](index);};

	inline float&		X_XYZ(const int Row, const int Col){return XYZ[0](Row, Col);};
	inline float&		Y_XYZ(const int Row, const int Col){return XYZ[1](Row, Col);};
	inline float&		Z_XYZ(const int Row, const int Col){return XYZ[2](Row, Col);};
	inline float&		X_XYZ(const long index){return XYZ[0](index);};
	inline float&		Y_XYZ(const long index){return XYZ[1](index);};
	inline float&		Z_XYZ(const long index){return XYZ[2](index);};

	inline float&		L_LAB(const int Row, const int Col){return LAB[0](Row, Col);};
	inline float&		A_LAB(const int Row, const int Col){return LAB[1](Row, Col);};
	inline float&		B_LAB(const int Row, const int Col){return LAB[2](Row, Col);};
	inline float&		L_LAB(const long index){return LAB[0](index);};
	inline float&		A_LAB(const long index){return LAB[1](index);};
	inline float&		B_LAB(const long index){return LAB[2](index);};
	//end by wtf

	inline int		gray(const long index);
	inline int		gray(const int Row, const int Col){return gray(cols()*Row + Col);};
			void	GrayImage(Raster<Tp> &image);


	inline long		rows() const {return _rows;};
	inline long		cols() const {return _cols;};
	inline long		size() const {return _rows*_cols;};

	inline bool		Valid(const int j, const int i)
						{if (j>=0&&j<rows()&&i>=0&&i<cols())return true;else return false;}

	inline long		Width(){return cols();};
	inline long		Height(){return rows();};

	virtual double	color_dis(const int index, const double mean_L, const double mean_U, const double mean_V);
	virtual void	Zero(void);

	// manipulation
	virtual void	convert_RGB_LUV(long selects=Lightness | Ustar | Vstar);			// converts RGB space to L*u*v space
	virtual void    convert_LUV_RGB(double *, int *, long select=Lightness | Ustar | Vstar);								// converts L*u*v space to RGB space
	
	//*******add by wtf
	virtual void    convert_RGB_YCrCb(void);
	virtual void    convert_YCrCb_RGB(double *pYCrCb, Tp *pRGB);
	virtual void    CreateYCrCbspace(void);
	virtual void    SetPixelInYCrCb(const int row, const int col,
									const int y, const int cr, const int cb);

	virtual void    convert_RGB_XYZ(void);
	virtual void    convert_XYZ_RGB(double *pXYZ, Tp *pRGB);
	virtual void    CreateXYZspace(void);
	virtual void    SetPixelInXYZ(const int row, const int col,
									const float x, const float y, const float z);

	virtual void    convert_RGB_LAB(bool scale=false);
	virtual void    convert_LAB_RGB(double *pLAB, Tp *pRGB, bool scale=false);
	virtual void    CreateLABspace(bool scale=false);
	virtual void    SetPixelInLAB(const int row, const int col,
									const float l, const float a, const float b, bool scale=false);

	virtual void    convert_XYZ_LAB(double *pXYZ, double *pLAB);
	virtual void    convert_LAB_XYZ(double *pLAB, double *pXYZ);

	virtual bool    IsColorImage();
	

	//*********end by wtf

	virtual void	histogram(long selects=Lightness | Ustar | Vstar);					// forms histogram of this color image
	virtual void	CreateLUVspace(void);
	virtual void	SetPixelInLUV(const int row, const int col,
								  const int l, const int u, const int v);
	virtual void	ToRGB(int nType=LUV2RGB);
	virtual void	Resize(Raster_Color<Tp>	&image_new, int ioption);

	// for color indexing
	virtual Matrix<int> Chromaticity(int r_g_b=RedColor);
	virtual void	MeanVariance(Matrix<int> &mask_label, irmm_Rect &bound,
								 const int num_label,
								 Matrix<double>	*pmeans,
								 Matrix<double> *pvariances,
								 Matrix<long> &counts);
	virtual void	Compute2DColor(Matrix<int> &mask_label, irmm_Rect &bound, const int num_label, Vector<double> *pvc_N);

protected:

	inline int		ShiftColorValues(const long index){return (((((int)R(index))<<8)+(int)G(index))<<8)+(int)B(index);};
	inline int		ShiftColorValues(const int Row, const int Col){return ShiftColorValues(rows() * cols() + Col);};


public:
	Raster<Tp>		RGB[3];
	Raster<int>		LUV[3];
	Raster<int>   YCrCb[3];       //add by wtf, for YCrCb color space
	Raster<float>   XYZ[3];
	Raster<float>   LAB[3];
	long			count;			//pixel count as row * col in this image
	Raster<uchar>	EdgeMap;		//edge map of this color image

protected:
	int			_rows, _cols;	//row number and column number of this image

	int          _n_colors;      //# of colors
	int**       _col_all;       //colors in the image (LUV)
	int*		_col0,*_col1,*_col2;      //to speed up
	int*		_col_RGB;       //colors in the image (RGB)
	int*        _col_index;     //color address
	int*        _col_misc;      //misc use in histogram and conversion
	int*        _m_colors;      //how many of each color
	int*        _col_remain;    //the map of remaining colors
	int*        _m_col_remain;  //table of remaining colors
	int         _n_col_remain;  //# of remaining colors
};


template<class Tp>
inline int Raster_Color<Tp>::gray(const long index)
{
	double a1,a2,a3;
	a1 = R(index);
	a2 = G(index);
	a3 = B(index);

	if (a1==a2 && a2==a3)
		return (int)R(index);
	else
		return (int)(0.299*a1+0.587*a2+0.114*a3);
};

template<class Tp1, class Tp2>
void CopyRasterColor(Raster_Color<Tp1>& x1,
						 Raster_Color<Tp2>& x2)
{
	x1.SetDimension(x2.rows(),x2.cols());
	for (int i=0; i<3; i++)
		CopyRaster(x1.RGB[i], x2.RGB[i]);
}

///////////////////////////////////////////////////////////
// operators
template<class Tp>
Raster_Color<Tp>& Raster_Color<Tp>::operator=(Raster_Color<uchar> &I)
{
	CopyRasterColor(*this,I);
	return *this;
}

template<class Tp>
Raster_Color<Tp>& Raster_Color<Tp>::operator=(Raster_Color<double> &I)
{
	CopyRasterColor(*this,I);
	return *this;
}

template<class Tp>
void Raster_Color<Tp>::Copy(Raster_Color<uchar> &I, irmm_Rect& rect)
{
	SetDimension(rect.Height(), rect.Width());
	for (int Row=rect.top; Row<=rect.bottom; Row++)
		for (int Col=rect.left; Col<=rect.right; Col++)
		{
			R(Row-rect.top,Col-rect.left) = (Tp)I.R(Row,Col);
			G(Row-rect.top,Col-rect.left) = (Tp)I.G(Row,Col);
			B(Row-rect.top,Col-rect.left) = (Tp)I.B(Row,Col);
		}
}

///////////////////////////////////////////////////////////
// inline functions
template<class Tp>
inline Tp& Raster_Color<Tp>::R(const int Row, const int Col)
{
	return RGB[0](Row, Col);
};

template<class Tp>
inline Tp& Raster_Color<Tp>::R(const long index)
{
	return RGB[0](index);
};

template<class Tp>
inline Tp& Raster_Color<Tp>::G(const int Row, const int Col)
{
	return RGB[1](Row, Col);
};

template<class Tp>
inline Tp& Raster_Color<Tp>::G(const long index)
{
	return RGB[1](index);
};

template<class Tp>
inline Tp& Raster_Color<Tp>::B(const int Row, const int Col)
{
	return RGB[2](Row, Col);
};

template<class Tp>
inline Tp& Raster_Color<Tp>::B(const long index)
{
	return RGB[2](index);
};

//==========================================================
// Class : Raster_Color
//==========================================================
//----------------------------------------------------------
// construction
//----------------------------------------------------------
////////////////////////////////////////////////////////////
template<class Tp>
Raster_Color<Tp>::Raster_Color()
{
	_rows = 0;
	_cols = 0;
	count = 0;
	_n_colors = 0;

	_col_all = NULL;
	_col0 = NULL;
	_col1 = NULL;
	_col2 = NULL;
	_col_index = NULL;
	_col_misc = NULL;
	_col_RGB = NULL;
	_m_colors = NULL;
}

////////////////////////////////////////////////////////////
template<class Tp>
Raster_Color<Tp>::~Raster_Color()
{
	int i = 0;
	Free_histogram_space();
}

template<class Tp>
void Raster_Color<Tp>::Free_histogram_space(void)
{
	if (_col_all != NULL)
	{
		delete [](_col_all[2]);
		delete [](_col_all[1]);
		delete [](_col_all[0]);
		delete [] _col_all;
	}
	if (_col_index != NULL)
		delete [] _col_index;
	if (_col_misc != NULL)
		delete [] _col_misc;
	if (_col_RGB != NULL)
		delete [] _col_RGB;
	if (_m_colors != NULL)
		delete [] _m_colors;
	_m_colors = NULL;
	_col_RGB = NULL;

	_col_misc = NULL;
	_col_index = NULL;
	_col_all = NULL;
}

template<class Tp>
void Raster_Color<Tp>::SetDimension(const int nrows, const int ncols)
{
	_rows = nrows;
	_cols = ncols;
	count = ((long)_rows) * _cols;
	_n_colors = 3;

	for (int i = 0; i < 3; i++)
		RGB[i].SetDimension(_rows, _cols);
}

template<class Tp>
void Raster_Color<Tp>::Zero(void)
{
	for (int i=0; i<3; i++)
	{
		LUV[i].Zero();
		RGB[i].Zero();
	}
}


//----------------------------------------------------------
// access
//----------------------------------------------------------
template<class Tp>
double Raster_Color<Tp>::color_dis(const int index, const double mean_L,
									 const double mean_U, const double mean_V)
{
	double fL, fU, fV;

	fL = L(index) - mean_L;
	fU = U(index) - mean_U;
	fV = V(index) - mean_V;
	return (fL * fL + fU * fU + fV * fV);
}

//---------------------------------------------------------------------
// manipulation

///////////////////////////////////////////////////////////////////////
// RGB to LUV conversion
// To gain speed the conversion works on a table of colors (_col_RGB[])
// rather than on the whole image
// Coefficient matrix for xyz and rgb spaces
static const int    RGB_LUV_XYZ[3][3] = { { 4125, 3576, 1804 },
                                  { 2125, 7154,  721 },
                                  {  193, 1192, 9502 } };
static const int	SEC_SIZE=64;  // 2^6

template<class Tp>
void Raster_Color<Tp>::convert_RGB_LUV(long selects)
{
	int		x, y, z, my_temp;
	double	l_star, u_star, v_star;
	double	u_prime, v_prime;
	 int temp_col, temp_index, temp_ind;
	 long j,k;

	int a00 = RGB_LUV_XYZ[0][0], a01 = RGB_LUV_XYZ[0][1], a02 = RGB_LUV_XYZ[0][2];
	int a10 = RGB_LUV_XYZ[1][0], a11 = RGB_LUV_XYZ[1][1], a12 = RGB_LUV_XYZ[1][2];
	int a20 = RGB_LUV_XYZ[2][0], a21 = RGB_LUV_XYZ[2][1], a22 = RGB_LUV_XYZ[2][2];

	int *A00 = new int[MAXV]; int *A01 = new int[MAXV]; int *A02 = new int[MAXV];
	int *A10 = new int[MAXV]; int *A11 = new int[MAXV]; int *A12 = new int[MAXV];
	int *A20 = new int[MAXV]; int *A21 = new int[MAXV]; int *A22 = new int[MAXV];

	for(j=0; j<MAXV; j++)
    {
		A00[j] = a00 * j; A01[j] = a01 * j; A02[j] = a02 * j;
		A10[j] = a10 * j; A11[j] = a11 * j; A12[j] = a12 * j;
		A20[j] = a20 * j; A21[j] = a21 * j; A22[j] = a22 * j;
    }

	double *my_pow = new double[MAXV];
	for (j = 0; j < MAXV; j++)
		my_pow[j]= (double)(116.0 * pow(j/255.0, 0.3333333) - 16);

	for (j = 0; j < _n_colors; j++)
    {
        temp_col = _col_RGB[j];
		int R = temp_col>>16;
		int G = (temp_col>>8) & 255;
		int B = temp_col & 255;

        x = A00[R] + A01[G] + A02[B];
        y = A10[R] + A11[G] + A12[B];
        z = A20[R] + A21[G] + A22[B];

        double  tval = (double)(y / 2550000.0); //Yn==1
		if ( tval >  Lt)
			l_star = my_pow[(int)(tval*255+0.5)];
        else
			l_star = 903.3 * tval;

        my_temp = x + 15 * y + 3 * z;
		if(my_temp)
		{
			u_prime = (double)(x << 2) / (double)(my_temp);
			v_prime = (double)(9 * y) / (double)(my_temp);
		}
		else
		{
			u_prime = 4.0;
			v_prime = 9.0 / 15.0;
		}

		tval = 13 * l_star;
        u_star = tval * (u_prime - Un_prime); // Un_prime = 0.1978
        v_star = tval * (v_prime - Vn_prime); // Vn_prime = 0.4683

		_col0[j] = (int)(l_star+0.5);
		if(u_star>0)
			_col1[j] = (int)(u_star+0.5);
        else
			_col1[j] = (int)(u_star-0.5);

		if(v_star>0)
			_col2[j] = (int)(v_star+0.5);
        else
			_col2[j] = (int)(v_star-0.5);
    }

	for (j = 0; j < count; j++)
    {
		temp_col = ShiftColorValues(j);
		temp_ind = _col_misc[temp_col>>6];
		for (k = temp_ind; k < temp_ind + SEC_SIZE; k++)
			if(_col_RGB[k] == temp_col)
			{
				temp_index = _col_index[j] = k;
				break;
			}
		L(j) = _col0[temp_index];
		U(j) = _col1[temp_index];
		V(j) = _col2[temp_index];
	}
	delete [] my_pow;
	delete [] A22; delete [] A21; delete [] A20;
	delete [] A12; delete [] A11; delete [] A10;
	delete [] A02; delete [] A01; delete [] A00;
	cerr<<":";
}

///////////////////////////////////////////////////////////////////
// LUV (final_T[]) to RGB (TI[]) conversion
template<class Tp>
void Raster_Color<Tp>::convert_LUV_RGB(double *pLUV, int *pRGB, long selects)
{
	int		i;
	double	R_G_B[3][3] = { {(double)3.2405, (double)-1.5371, (double)-0.4985 },
                            {(double)-0.9693, (double)1.8760, (double)0.0416 },
                            {(double)0.0556, (double)-0.2040, (double)1.0573 }};

	// this condition is always true
	if (selects & Lightness && selects & Ustar && selects & Vstar)
	{
		if(pLUV[0] < 0.1)
        {
			pRGB[0]=0;
			pRGB[1]=0;
			pRGB[2]=0;
		}
		else
		{
			double my_x, my_y, my_z;
			if (pLUV[0] < 8.0)
				my_y = Yn * pLUV[0] / 903.3;
			else
				my_y = Yn * pow((pLUV[0] + 16.0) / 116.0, 3);

			double u_prime = pLUV[1] / (13 * pLUV[0]) + Un_prime;
			double v_prime = pLUV[2] / (13 * pLUV[0]) + Vn_prime;

			my_x = 9 * u_prime * my_y / (4 * v_prime);
			my_z = (12 - 3 * u_prime - 20 * v_prime) * my_y / (4 * v_prime);

			pRGB[0] = int((R_G_B[0][0]*my_x + R_G_B[0][1]*my_y + R_G_B[0][2]*my_z)*255.0);
			pRGB[1] = int((R_G_B[1][0]*my_x + R_G_B[1][1]*my_y + R_G_B[1][2]*my_z)*255.0);
			pRGB[2] = int((R_G_B[2][0]*my_x + R_G_B[2][1]*my_y + R_G_B[2][2]*my_z)*255.0);

			for (i = 0; i < 3; i++)
			{
				if (pRGB[i] > 255) pRGB[i] = 255;
				else if (pRGB[i] < 0) pRGB[i] = 0;
			}
		}
    }
	else
    {
		pRGB[0] = (int)pLUV[0];
		pRGB[1] = (int)pLUV[1];
		pRGB[2] = (int)pLUV[2];
	}
}

template<class Tp>
void Raster_Color<Tp>::ToRGB(int nType)
{
	int i,j,pRGB[3];
	double pLUV_YCrCb[3];

	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			if ( nType==LUV2RGB )
			{
				pLUV_YCrCb[0] = (double)L(j,i);
				pLUV_YCrCb[1] = (double)U(j,i);
				pLUV_YCrCb[2] = (double)V(j,i);
					
				convert_LUV_RGB(pLUV_YCrCb, pRGB);
			}
			else
			{
				pLUV_YCrCb[0] = (double)Y(j,i);
				pLUV_YCrCb[1] = (double)Cr(j,i);
				pLUV_YCrCb[2] = (double)Cb(j,i);

				convert_YCrCb_RGB(pLUV_YCrCb, (Tp *)pRGB);
			}			

			R(j,i) = (Tp)(pRGB[0]);
			G(j,i) = (Tp)(pRGB[1]);
			B(j,i) = (Tp)(pRGB[2]);
		}
}


//////////////////////////////////////////////////////////////
// 3-D Histogram computation
// Implement a trade-off between speed and required memory
// Used in 3-D histogram computation
static const int FIRST_SIZE=262144; // 2^18

template<class Tp>
void Raster_Color<Tp>::histogram(long selects)
{
	// free assigned space
	Free_histogram_space();

	int *first_tab= new int[FIRST_SIZE];
	int **third_tab;
	int *fourth_tab;
	int *fifth_tab=new int[SEC_SIZE];
	uchar oR,oG,oB;

	_n_colors=0;

	long k,j,p,r;
	int temp_ind, sec_ind, third_ind;

	int first_contor=0, third_contor=0;

	memset(first_tab, 0, sizeof(int) * FIRST_SIZE);
	_col_misc= new int[FIRST_SIZE];
	memset(_col_misc, 0, sizeof(int) * FIRST_SIZE);

	//first_tab -> how many
	for(k = 0; k < count; k++)
    {
		oR = (uchar)R(k);
		oG = (uchar)G(k);
		oB = (uchar)B(k);

		temp_ind = (((oR<<8)+oG)<<2)+(oB>>6);
		first_tab[temp_ind]++;
    }
	//_col_misc -> memo position
	for(k = 0; k < FIRST_SIZE; k++)
		if(first_tab[k])
		{
			_col_misc[k]=first_contor;
			first_contor++;
		}
	//contors
	fourth_tab = new int[first_contor];
	memset(fourth_tab, 0, sizeof(int) * first_contor);
	//tab of pointers to reduced colors
	third_tab = new int *[first_contor];
	first_contor = 0;
	for (k =0 ; k < FIRST_SIZE; k++)
		if(first_tab[k])
		{
			third_tab[first_contor] = new int[first_tab[k]];
			first_contor++;
		}

	for (k = 0; k < count; k++)
    {
		oR = (uchar)R(k);
		oG = (uchar)G(k);
		oB = (uchar)B(k);

		temp_ind = (((oR<<8)+oG)<<2)+(oB>>6);
		sec_ind = oB & 63;
		third_ind=_col_misc[temp_ind];
		third_tab[third_ind][fourth_tab[third_ind]]  = sec_ind;
		fourth_tab[third_ind]++;
	}

	for (k = 0; k < first_contor; k++)
    {
		memset(fifth_tab, 0, sizeof(int)*SEC_SIZE);
		for (j=0;j<fourth_tab[k];j++)
			fifth_tab[third_tab[k][j]]++;
		for (j = 0; j < SEC_SIZE; j++)
			if(fifth_tab[j])
				_n_colors++;
    }

	_col_RGB=new int[_n_colors];
	_m_colors=new int[_n_colors];

	k = 0; p = 0;
	for (r = 0; r < FIRST_SIZE; r++)
		if(first_tab[r])
		{
			memset(fifth_tab,0,sizeof(int)*SEC_SIZE);
			for(j = 0; j < fourth_tab[k]; j++)
				fifth_tab[third_tab[k][j]]++;
			_col_misc[r] = p;
			for(j = 0; j < SEC_SIZE; j++)
				if(fifth_tab[j])
				{
					_col_RGB[p] = (r<<6) + j;
					_m_colors[p] = fifth_tab[j];
					p++;
				}
			delete [] third_tab[k];
			k++;
		}

	delete [] third_tab;
	delete [] fourth_tab;
	delete [] fifth_tab;
	delete [] first_tab;

	_col_all = new int*[3];
	_col0=_col_all[0] = new int[_n_colors];
	_col1=_col_all[1] = new int[_n_colors];
	_col2=_col_all[2] = new int[_n_colors];
	_col_index = new int[count];
	cerr<<":";
}

template<class Tp>
void Raster_Color<Tp>::CreateLUVspace(void)
{
	long	i;

	// get histogram of this image
	histogram();
	for (i = 0; i < 3; i++)
		LUV[i].SetDimension(Height(), Width());
	convert_RGB_LUV();
}

// for color indexing
template<class Tp>
Matrix<int> Raster_Color<Tp>::Chromaticity(int r_g_b)
{
	double		d_chr;
	Matrix<int>	chr(rows(),cols());
	int			i,j,R_G_B;

	for (j=0; j<rows(); j++)
		for (i=0; i<cols(); i++)
		{
			R_G_B = (int)(R(j,i)+G(j,i)+B(j,i));
			if (R_G_B == 0)
				d_chr = 1.0;
			else
			{
				if (r_g_b == RedColor)
					d_chr = ((double)R(j,i))/((double)(R_G_B));
				else
					d_chr = ((double)G(j,i))/((double)(R_G_B));
			}

			chr(j,i) = (int)(d_chr*255);
		}
	return chr;
}

template<class Tp>
void Raster_Color<Tp>::SetPixelInLUV(const int row, const int col,
							const int l, const int u, const int v)
{
	double	luv[3];
	int		rgb[3];

	luv[0] = double(l);
	luv[1] = double(u);
	luv[2] = double(v);

	convert_LUV_RGB(luv, rgb);
	R(row,col) = (Tp)rgb[0];
	G(row,col) = (Tp)rgb[1];
	B(row,col) = (Tp)rgb[2];
}

template<class Tp>
void Raster_Color<Tp>::MeanVariance(Matrix<int> &mask_label, irmm_Rect &bound,
								 const int num_label,
								 Matrix<double>	*pmeans,
								 Matrix<double> *pvariances,
								 Matrix<long> &counts)
{
	int i,j,m,n,k;

	if (pmeans == NULL) return;
	long *pnum_pixel=new long[num_label];

	for (i=0; i<num_label; i++)
	{
		pnum_pixel[i] = 0;
		pmeans[i].SetDimension(3,1);
		pmeans[i].InitValue(0.0);
		if (pvariances != NULL)
		{
			pvariances[i].SetDimension(3,3);
			pvariances[i].InitValue(0.0);
		}
	}
	counts.SetDimension(num_label,1);

	// compute the mean value
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if (mask_label(j,i)>0 && mask_label(j,i)<=num_label)
			{
				k = mask_label(j,i)-1;
				pnum_pixel[k]++;
				pmeans[k](0) += (double)L(j,i);
				pmeans[k](1) += (double)U(j,i);
				pmeans[k](2) += (double)V(j,i);
			}
		}
	for (i=0; i<num_label; i++)
	{
		if (pnum_pixel[i]>0)
		{
			pmeans[i](0) = pmeans[i](0)/((double)pnum_pixel[i]);
			pmeans[i](1) = pmeans[i](1)/((double)pnum_pixel[i]);
			pmeans[i](2) = pmeans[i](2)/((double)pnum_pixel[i]);
		}
		else
		{
			pmeans[i](0) = 0.0;
			pmeans[i](1) = 0.0;
			pmeans[i](2) = 0.0;
		}
	}

	if (pvariances != NULL)
	{
		for (j=bound.top; j<=bound.bottom; j++)
			for (i=bound.left; i<=bound.right; i++)
			{
				if (mask_label(j,i)>0 && mask_label(j,i)<=num_label)
				{
					k = mask_label(j,i)-1;
					for (n=0; n<3; n++)
						for (m=0; m<3; m++)
							pvariances[k](n,m) += ((double)(LUV[n](j,i)-pmeans[k](n)))*
												  ((double)(LUV[m](j,i)-pmeans[k](m)));
				}
			}
		for (i=0; i<num_label; i++)
		{
			for (n=0; n<3; n++)
				for (m=0; m<3; m++)
				{
					if (pnum_pixel[i]>1)
						pvariances[i](n,m) = pvariances[i](n,m)/(pnum_pixel[i]-1);
				}
			counts(i) = (long)pnum_pixel[i];
		}
		// make sure the determinate of covariance matrix is not zero
		for (i=0; i<num_label; i++)
		{
			if (pvariances[i].Det()==0)
				pvariances[i].Identity();
		}
	}
	delete []pnum_pixel;
}

template<class Tp>
void Raster_Color<Tp>::Compute2DColor(Matrix<int> &mask_label, irmm_Rect &bound,const int num_label, Vector<double> *pvc_N)
{
	int				i,j,m,k;
	double			din;
	Vector<double>	vc_weights,vc_temp(3);

	if (pvc_N == NULL) return;

	for (i=0; i<num_label; i++)
	{
		pvc_N[i].SetRowVector(3);
		pvc_N[i].InitValue(0.0);
	}
	vc_weights.SetRowVector(num_label);
	vc_weights.InitValue(0.0);

	// compute the major direction
	for (j=bound.top; j<=bound.bottom; j++)
		for (i=bound.left; i<=bound.right; i++)
		{
			if (mask_label(j,i)>=0 && mask_label(j,i)<num_label)
			{
				k = mask_label(j,i);

				vc_temp(0) = (double)R(j,i);
				vc_temp(1) = (double)G(j,i);
				vc_temp(2) = (double)B(j,i);

				// get intensity
				din = vc_temp.Norm(2.0);
				// normalize it
				if (din>0.0)
				{
					for (m=0; m<3; m++)
					{
						vc_temp(m) = vc_temp(m)/din;
						pvc_N[k](m) = (pvc_N[k](m)*vc_weights(k)+din*vc_temp(m))/(vc_weights(k)+din);
					}
					vc_weights(k) += din;
				}
			}
		}

	for (k=0; k<num_label; k++)
	{
		din = pvc_N[k].Norm(2.0);
		if (din>0.0)
			pvc_N[k] = pvc_N[k]/din;
	}
}


// create a new image with the size specified in image_new
//  option 0: do a simple sampling, no
template<class Tp>
void Raster_Color<Tp>::Resize(Raster_Color<Tp>	&image_new, int ioption)
{
	double	delta_row,delta_col;
	int		i,j,row,col;

	if (image_new.rows()*cols()<=1 || rows()*cols()<=1)
		return;

	if (ioption == 0)
	{
		delta_row = ((double)(rows()-1))/(image_new.rows()-1);
		delta_col = ((double)(cols()-1))/(image_new.cols()-1);

		for (j=0; j<image_new.rows(); j++)
			for (i=0; i<image_new.cols(); i++)
			{
				row = IntClose(delta_row*j);
				col = IntClose(delta_col*i);

				image_new.R(j, i) = R(row, col);
				image_new.G(j, i) = G(row, col);
				image_new.B(j, i) = B(row, col);
			}
	}
}


template<class Tp>
void Raster_Color<Tp>::LoadPpmPGm(const char *filename)
{
	fstream fin(filename, ios::in | ios::binary | ios::_Nocreate );
	if( fin.fail())
	{
		cerr << "Fail to open input file " << filename << " in ColorImage::read()" << endl;
		return;
	}

	char	magicNumber[2],comment[281];
	int		MN,xmax,ymax,image_max;

	fin.getline(magicNumber, 10);
	MN=atoi(magicNumber+1);
	while(fin.peek() == '#' )  fin.getline( comment, 81 );  //remove comments
	while( !IsDigit( (int)fin.peek() ) )
		fin.get(magicNumber[1]);
	fin >> xmax >> ymax >> image_max;
	if(xmax!=cols() || ymax!=rows())
		SetDimension(ymax, xmax);
	while(fin.peek()!='\n') fin.get();
		fin.get(); //remove extra newline

	unsigned char * tempbuffer;
	long			i;

	switch( MN)
	{
	case iffPpmAscii:	//ASCII PPM
						int x,y,z;
						for (i=0; i<size(); i++)
						{
							fin >> x >> y >> z;
							R(i) = (Tp)x;
							G(i) = (Tp)y;
							B(i) = (Tp)z;
						}
						break;
	case iffPpmBinary:	//BINARY PPM
						tempbuffer = new unsigned char[size()*3];
						fin.read((char *)tempbuffer, size()*3);
			for (i=0; i<size(); i++)
			{
				R(i) = (Tp)tempbuffer[i*3];
				G(i) = (Tp)tempbuffer[i*3+1];
				B(i) = (Tp)tempbuffer[i*3+2];
			}
			delete [] tempbuffer;
			break;
	case iffPgmAscii: //ASCII PGM
			int input;
			for (i=0; i<size(); i++)
			{
				fin >> input;
				R(i) = (Tp)input;
				G(i) = R(i);
				B(i) = R(i);
			}
			break;
	case iffPgmBinary: //BINARY PGM
			tempbuffer = new unsigned char[size()];
			fin.read((char *)tempbuffer, size());
			for (i=0; i<size(); i++)
			{
				R(i) = (Tp)tempbuffer[i];
				G(i) = R(i);
				B(i) = R(i);
			}
			delete []tempbuffer;
			break;
	}
	fin.close();
}

template<class Tp>
bool Raster_Color<Tp>::SavePpmPGm(const char * filename, int iff)
{
	int i;

	assert(size()>0);
	fstream fout( filename, ios::out | ios::binary );
	if( fout.fail() )
	{
		cerr << "Error opening output file ColorImage::Write()" << endl;
		return false;
	}
	fout << "P" << iff << endl << cols() << ' ' << rows() << endl << 255 << endl;
	switch(iff)
	{
	case iffPpmAscii: //ASCII PPM
			for(i=0; i<size(); i++ )
				fout << int(R(i)) << ' ' << int(G(i))
					 << ' ' << int(B(i)) <<' ';
			break;
	case iffPpmBinary: //BINARY PPM
			char* tempbuffer;
			tempbuffer = new char[size()*3];
			for(i=0; i<size(); i++)
			{
				tempbuffer[i*3]		= (char)R(i);
				tempbuffer[i*3+1]	= (char)G(i);
				tempbuffer[i*3+2]	= (char)B(i);
			}
			fout.write(tempbuffer,(size()*3));
			delete []tempbuffer;
			break;
	case iffPgmAscii: //ASCII PGM
			for(i=0; i<size(); i++)
				fout << (int) gray(i) << ' ';
			break;
	case iffPgmBinary: //BINARY PGM
			tempbuffer = new char[size()];
			for(i=0; i<size(); i++)
				tempbuffer[i] = (char)gray(i);
			fout.write(tempbuffer, size());
			delete []tempbuffer;
			break;
	}
	return true;
}

template<class Tp>
void Raster_Color<Tp>::GrayImage(Raster<Tp> &image)
{
	int	i,j;

	image.SetDimension(rows(), cols());
	for (j=0; j<_rows; j++)
		for (i=0; i<_cols; i++)
			image(j, i) = gray(j, i);
}

//add by wtf for integrating PS code
template<class Tp>
void Raster_Color<Tp>::convert_RGB_YCrCb()
{
	// RGB => YCrCb
	    // Y=0.299*R + 0.587*G + 0.114*B
		// Cr=(R-Y)*0.713 + 128
		// Cb=(B-Y)*0.564 + 128

	double Y=0.0, Cr=0.0, Cb=0.0;

	for ( int i=0; i<_rows; i++ )
		for ( int j=0; j<_cols; j++ )
			{
				Y = RGB[0](i,j) * 0.299 + RGB[1](i,j) * 0.587 + RGB[2](i,j) * 0.114;
				Cr = ((double)RGB[0](i,j) - Y) * 0.713 + 128.0;
				Cb = ((double)RGB[2](i,j) - Y) * 0.564 + 128.0;

				if ( Y<0 )
				{
					YCrCb[0](i,j) = 0;
				}
				else if ( Y>255 )
				{
					YCrCb[0](i,j) = 255;
				}
				else
				{
					YCrCb[0](i,j) = (uchar)Y;
				}

				if ( Cr<0 )
				{
					YCrCb[1](i,j) = 0;
				}
				else if ( Cr>255 )
				{
					YCrCb[1](i,j) = 255;
				}
				else
				{
					YCrCb[1](i,j) = (uchar)Cr;
				}

				if ( Cb<0 )
				{
					YCrCb[2](i,j) = 0;
				}
				else if ( Cb>255 )
				{
					YCrCb[2](i,j) = 255;
				}
				else
				{
					YCrCb[2](i,j) = (uchar)Cb;
				}

              }
	
}

template<class Tp>
void Raster_Color<Tp>::convert_YCrCb_RGB(double *pYCrCb, Tp *pRGB)
{
	//YCrCb => RGB
	    // R=Y + 1.403*(Cr - 128)
		// G=Y - 0.344*(Cr - 128) - 0.714*(Cb - 128)
		// B=Y + 1.773*(Cb - 128)
	double R=0.0, G=0.0, B=0.0;

	R = pYCrCb[0] + (pYCrCb[1] - 128) * 1.403;
	G = pYCrCb[0] - (pYCrCb[1] - 128) * 0.344 - (pYCrCb[2] - 128) * 0.714;
	B = pYCrCb[0] + (pYCrCb[2] - 128) * 1.773;

	if ( R<0 )
	{
		pRGB[0] = 0;
	}
	else if ( R>255 )
	{
		pRGB[0] = 255;
	}
	else
	{
		pRGB[0] = (int)R;
	}

	if ( G<0 )
	{
		pRGB[1] = 0;
	}
	else if ( G>255 )
	{
		pRGB[1] = 255;
	}
	else
	{
		pRGB[1] = (int)G;
	}

	if ( B<0 )
	{
		pRGB[2] = 0;
	}
	else if ( B>255 )
	{
		pRGB[2] = 255;
	}
	else
	{
		pRGB[2] = (int)B;
	}
}

template<class Tp>
void Raster_Color<Tp>::CreateYCrCbspace()
{	
	// get histogram of this image
	histogram();
	for (int i = 0; i < 3; i++)
		YCrCb[i].SetDimension(Height(), Width());

	convert_RGB_YCrCb();

}

template<class Tp>
void Raster_Color<Tp>::SetPixelInYCrCb(const int row, const int col, const int y, const int cr, const int cb)
{
	double	ycrcb[3];
	Tp		rgb[3];

	ycrcb[0] = double(y);
	ycrcb[1] = double(cr);
	ycrcb[2] = double(cb);

	convert_YCrCb_RGB(ycrcb, rgb);

	R(row,col) = (Tp)rgb[0];
	G(row,col) = (Tp)rgb[1];
	B(row,col) = (Tp)rgb[2];
}


template<class Tp>
void Raster_Color<Tp>::convert_RGB_XYZ()
{
	
	double r=0.0, g=0.0, b=0.0;
	double r_prime=0.0, g_prime=0.0, b_prime=0.0;
	double X=0.0, Y=0.0, Z=0.0;
	const double power = 2.4;


	
	for ( int i=0; i<_rows; i++ )
		for ( int j=0; j<_cols; j++ )
		{
			r_prime = RGB[0](i, j) / 255.;
			g_prime = RGB[1](i, j) / 255.;
			b_prime = RGB[2](i, j) / 255.;

			r = ( r_prime <= 0.04045 ) ? r = r_prime / 12.92 : pow((double)(( r_prime + 0.055 ) / 1.055), power);

			g = ( g_prime <= 0.04045 ) ? g_prime / 12.92 : pow((double)(( g_prime + 0.055 ) / 1.055), power);

			b = ( b_prime <= 0.04045 ) ? b_prime / 12.92 : pow((double)(( b_prime + 0.055 ) / 1.055), power);

			r *= 100;
			g *= 100;
			b *= 100;
			X = 0.412424 * r + 0.357579 * g + 0.180464 * b;
			Y = 0.212656 * r + 0.715158 * g + 0.0721856 * b;
			Z = 0.0193324 * r + 0.119193 * g + 0.950444 * b;

			XYZ[0](i,j) = static_cast<float>(((X * 10000) + .5)/10000.);
			XYZ[1](i,j) = static_cast<float>(((Y * 10000) + .5)/10000.);
			XYZ[2](i,j) = static_cast<float>(((Z * 10000) + .5)/10000.);

		}
}

template<class Tp>
void Raster_Color<Tp>::convert_XYZ_RGB(double *pXYZ, Tp *pRGB)
{
	double r, g, b;
	double R, G, B;

	r = 3.24071 * pXYZ[0] + (-1.53726) * pXYZ[1] + (-0.498571) * pXYZ[2];
	g = (-0.969258) * pXYZ[0] + 1.87599 * pXYZ[1] +	0.0415557 * pXYZ[2];
	b = 0.0556352 * pXYZ[0] +	(-0.203996) * pXYZ[1] + 1.05707 * pXYZ[2];

	r = static_cast<int>(((r * 100) + .5)/10000.);
	g = static_cast<int>(((g * 100) + .5)/10000.);
	b = static_cast<int>(((b * 100) + .5)/10000.);

	// 
	r = ( r < 0 ) ? 0. : r;
	r = ( r > 1 ) ? 1. : r;
	g = ( g < 0 ) ? 0. : g;
	g = ( g > 1 ) ? 1. : g;
	b = ( b < 0 ) ? 0. : b;
	b = ( b > 1 ) ? 1. : b;

	if ( r <= 0.0031308 )
		R = 12.92 * r;
	else
		R = 1.055 * pow((double)r, (double)(1.0/2.4)) - 0.055;

	if ( g <= 0.0031308 )
		G = 12.92 * g;
	else
		G = 1.055 * pow((double)g, (double)1.0/2.4) - 0.055;

	if ( b <= 0.0031308 )
		B = 12.92 * b;
	else
		B = 1.055 * pow((double)b, (double)(1.0/2.4)) - 0.055;

	pRGB[0] = static_cast<int> (floor(R * 255+.5));
	pRGB[1] = static_cast<int> (floor(G * 255+.5));
	pRGB[2] = static_cast<int> (floor(B * 255+.5));

	pRGB[0] = ( pRGB[0] > 255 ) ? 255 : pRGB[0];
	pRGB[1] = ( pRGB[1] > 255 ) ? 255 : pRGB[1];
	pRGB[2] = ( pRGB[2] > 255 ) ? 255 : pRGB[2];
}

template<class Tp>
void Raster_Color<Tp>::CreateXYZspace()
{
	histogram();
	for (int i = 0; i < 3; i++)
		XYZ[i].SetDimension(Height(), Width());

	convert_RGB_XYZ();
}

template<class Tp>
void Raster_Color<Tp>::SetPixelInXYZ(const int row, const int col, const float x, const float y, const float z)
{
	double	xyz[3];
	Tp		rgb[3];

	xyz[0] = double(x);
	xyz[1] = double(y);
	xyz[2] = double(z);

	convert_XYZ_RGB(xyz, rgb);

	R(row,col) = (Tp)rgb[0];
	G(row,col) = (Tp)rgb[1];
	B(row,col) = (Tp)rgb[2];
}

template<class Tp>
void Raster_Color<Tp>::convert_XYZ_LAB(double *pXYZ, double *pLAB)
{
	double Xr, Yr, Zr, x, y, z;
	double eps, k, fx, fy, fz;
	const bool intentCIE = false;

	eps = intentCIE ? 216./24389 : 0.008856;
	k = intentCIE ? 24389./27 : 903.3;
	// CIE D65 Reference White
	Xr = 95.047;
	Yr = 100.00;
	Zr = 108.883;

	x = pXYZ[0] / Xr;
	y = pXYZ[1] / Yr;
	z = pXYZ[2] / Zr;

	fx = x > eps ? pow((double)x, (double)(1./3.)) : (k*x+16)/116.;
	fy = y > eps ? pow((double)y, (double)(1./3.)) : (k*y+16)/116.;
	fz = z > eps ? pow((double)z, (double)(1./3.)) : (k*z+16)/116.;

	pLAB[0] = 116*fy - 16;
	pLAB[1] = 500 * (fx-fy);
	pLAB[2] = 200 * (fy-fz);
}

template<class Tp>
void Raster_Color<Tp>::convert_LAB_XYZ(double *pLAB, double *pXYZ)
{
	double Xr, Yr, Zr, x, y, z;
	double eps, k, fx, fy, fz;
	const bool intentCIE = false;

	eps = intentCIE ? 216./24389 : 0.008856;
	k = intentCIE ? 24389./27 : 903.3;

	Xr = 95.047;
	Yr = 100.00;
	Zr = 108.883;


	y= ( pLAB[0] > k*eps ) ? pow((double)((pLAB[0]+16)/116), 3.) : pLAB[0]/k;
	fy = (y > eps) ? (pLAB[0]+16)/116 : (k*y+16)/116;
	fz = fy - pLAB[2]/200;
	fx = pLAB[1] / 500 + fy;
	x = ( pow((double)fx, 3.) > eps ) ? pow((double)fx, 3.) : (116 * fx - 16) / k;
	z = ( pow((double)fz, 3.) > eps ) ? pow((double)fz, 3.) : (116 * fz - 16) / k;

	pXYZ[0] = Xr * x;
	pXYZ[1] = Yr * y;
	pXYZ[2] = Zr * z;
}

template<class Tp>
void Raster_Color<Tp>::convert_RGB_LAB(bool scale)
{
	convert_RGB_XYZ();

	double xyz[3], lab[3];
	
	for ( int i=0; i<_rows; i++ )
		for ( int j=0; j<_cols; j++ )
		{
			xyz[0] = XYZ[0](i,j);
			xyz[1] = XYZ[1](i,j);
			xyz[2] = XYZ[2](i,j);

			convert_XYZ_LAB(xyz, lab);

			if ( scale )
			{
				lab[0] /= 100.;
				lab[1] /= 100.;
				lab[2] /= 110.;
			}

			LAB[0](i, j) = (float)(lab[0]); 
			LAB[1](i, j) = (float)(lab[1]); 
			LAB[2](i, j) = (float)(lab[2]); 			
		}

}

template<class Tp>
void Raster_Color<Tp>::convert_LAB_RGB(double *pLAB, Tp *pRGB,bool scale)
{
	if ( scale )
	{
		pLAB[0] *= 100.;
		pLAB[1] *= 100.;
		pLAB[2] *= 110.;
	}

	double xyz[3];
	convert_LAB_XYZ(pLAB, xyz);
	convert_XYZ_RGB(xyz, pRGB);

}

template<class Tp>
void Raster_Color<Tp>::CreateLABspace(bool scale)
{
	histogram();
	for (int i = 0; i < 3; i++)
	{
		XYZ[i].SetDimension(Height(), Width());
		LAB[i].SetDimension(Height(), Width());
	}

	convert_RGB_LAB(scale);
}

template<class Tp>
void Raster_Color<Tp>::SetPixelInLAB(const int row, const int col, const float l, const float a, const float b, bool scale )
{
	Tp rgb[3];
	double lab[3];

	lab[0] = l;
	lab[1] = a;
	lab[2] = b;

	convert_LAB_RGB( lab, rgb, scale);

	R(row,col) = (Tp)rgb[0];
	G(row,col) = (Tp)rgb[1];
	B(row,col) = (Tp)rgb[2];

}

template<class Tp>
bool Raster_Color<Tp>::IsColorImage()
{

	unsigned int	i, j;
	Tp			byR,byG,byB;
	for (j=0; j<_rows; j++)
		for (i=0; i<_cols; i++)
		{
			byR = RGB[0](j, i);
			byG = RGB[1](j, i);
			byB = RGB[2](j, i);

			if (abs(byR-byG)>5 || abs(byR-byB)>5 ||abs(byG-byB)>5)
				return true;	// any three components in color are different
		}

	return false;
}

//end by wtf

#endif
